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/""N ' In this paper we consider a resonant injection- locked frequency divider which is of interest in 

electronics, and we investigate the frequency locking phenomenon when varying the amplitude 
i-C • and frequency of the injected signal. We study both analytically and numerically the structure 

of the Arnol'd tongues in the frequency-amplitude plane. In particular, we provide exact 
analytical formulae for the widths of the tongues, which correspond to the plateaux of the 
devil's staircase picture. The results account for numerical and experimental findings presented 
in the literature for special driving terms and, additionally, extend the analysis to a more 
general setting. 
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Abstract 



The locking of oscillators onto subharmonics of the driving frequency (known as frequency locking 
or frequency demultiplication) has been well known in electronics since the work of van der Pol 
and van der Mark [22]; see also [12]. In the frequency-amplitude plane, locking occurs in distorted 
wedge-shaped regions (Arnol'd tongues) with apices corresponding to rational values on the (scaled) 
frequency axis. If one plots the ratio of the driving frequency to the output frequency versus the 
driving frequency, one obtains a so-called devil's staircase, i.e. a self-similar fractal object, where 
the qualitative structure is replicated at higher levels of resolution, with plateaux corresponding 
to rational values of the ratio. 

The frequency locking phenomenon, the existence of the Arnold tongues, and the devil's stair- 
case structure have been proved rigorously in some mathematical models, such as the circle map [3], 
and studied numerically or experimentally for several electronic circuits, such as the van der Pol 
equation [9, 17], the Josephson junction [1, 13, 20], the Chua circuit [18] among others. 

In this paper we are interested in studying both analytically and numerically an electronic 
circuit, namely a resonant injection-locked frequency divider (ILFD), first considered in [15]. In [16] 
a differential equation was introduced to describe the circuit and it was shown that the numerical 
integration of the equation reliably reproduces the experimental data. 

In [4] , the differential equation describing the ILFD was studied with the purpose of explaining 
analytically the appearance of frequency locking. In particular, the full differential equation in 
question was shown to be of the form u" + u'h(u) + k(u) + fi^!(u, u' , t) — 0, where h(u) and k(u) 
are even and odd functions of u, respectively, the prime denotes differentiation with respect to t, 
and fi is the amplitude of the perturbation 'J, which is taken to be periodic in t, with frequency lo, 
i.e. ^f(u,u',t) = ^(u, u', t + 2h/lo). If the drive is purely sinusoidal, as in [16], the Fourier series 
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expansion of ^>(u, u' , t), 

contains only the first harmonics v = ±1 (i.e. ty v (u,u') — for \v\ > 1). 

When p = 0, the system is unperturbed, and the differential equation is of a particular form 
known as the Lienard equation [10, 9]; the best known example of this type is the van der Pol 
equation. Under suitable assumptions on h and k, the Lienard equation admits a globally-attracting 
limit cycle. 

The phenomenon of frequency locking manifests itself in the ILFD when the ratio of the driv- 
ing frequency u> to the output frequency Q is plotted against lu. When lu is close to a rational 
multiple p = p/q of the frequency f^o of the limit cycle of the unperturbed system, then f2 is 
fixed such that u> = pVl. Therefore the plot has a devil's staircase structure [15], with plateaux 
corresponding to rational values of the ratio lj/Q. If lj/Q — p/q one says that there is a resonance 
(or synchronisation) of order p:q. For purely sinusoidal perturbations \f, such as those considered 
explicitly in [15, 16], the main plateaux correspond to even values of p (a physical argument was 
given in [23]). The perturbation theory approach taken in [4] successfully explains the experimental 
observations, by computing quantitatively the way in which the widths of the plateaux depend on 
the amplitude of the perturbation p, assumed small. 

In an alternative visualisation of the situation, a two-dimensional plot showing where locking 
takes place is constructed in the (to,p) plane. The Arnold tongues have widths and centre-lines 
which vary as some (explicitly computable) integer power of p [4]. The experimentally-observed 
dominance of tongues for which the ratio lu/^Iq is close to an even integer can be explained by the 
fact that only these tongues have a width of order p: all other tongues grow in width as some higher 
power of p. Specifically, if p G Q and Au(p) — {lu : u>/Q = p} is the width of the corresponding 
locking interval at fixed p, it was proved that, for sinusoidal perturbations, 

Aw(2n/fc) = 0(p k ), Auj{{2n-l)/k) = 0(p 2k ) (1.1) 

for all k,n G IN such that 2n/k and (2n + l)/k are irreducible fractions. The centre-lines are 
vertical for p = 2n and bend away from the vertical by a distance of order p 2 for all other values 
of p. In [4] we also stressed that the property (1.1) strongly depends on the particular form of the 
drive, more precisely on the fact that, as in [15, 16], the driving was taken to contain only the first 
harmonics. 

More generally, one can express Aiv(p) as a power series (perturbation series) in p, 

oo 

Au;(p) = J2^^Mp). (1.2) 

k=i 

If ko € IN is the first integer such that Auuk Q {p) ^ then from (1.2) one obtains Auj(p) = 
p k "A k Mp) + O(p k0+1 ). 

The convergence of the perturbation series for p small enough — yielding analyticity in p in a 
neighbourhood of the origin — was discussed and proved in [4]; see also [8]. Hence, by keeping only 
the lowest order terms means that in (1.1) we are looking at the leading contributions, without 
making any uncontrolled truncation. The coefficients AkLu(p) are given in the form of suitable 
integrals. However, it is not possible to reduce this computation to the integration of elementary 
functions, because the integrands involve functions which are known only numerically. Thus, the 
computation of the integrals requires some work, which we also discuss in this paper. 

A first order analysis of the locking intervals (in the same spirit as in [4]) is also performed 
in [2], where only sinusoidal perturbations are considered; in particular no prediction is made 
for resonances p : q with p/q ^ 21N, as this would require a higher order analysis. More general 
perturbations are considered in [14], where a different approach is followed. However, this involves 
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approximations which, ultimately, correspond to a first order analysis. By contrast, the analysis 
performed in [8] and further developed here allows us to go to arbitrary perturbation orders, with a 
control on the remainder. Thus, not only one can find an exact analytical expression for the leading 
order of the locking interval of any resonance, but in principle one can also compute any locking 
interval within any desired accuracy. In [7], the lack of accurate analytical methods to predict 
the locking range was deplored: in our opinion our analysis, which makes no approximation, fills 
this gap. Of course, for practical purposes, the computation of the locking interval for any given 
resonance requires solving numerically some integrals (which become increasingly complicated as 
the perturbation order increases). It would be desirable to have a formula for the locking interval 
in terms of the parameters a and (3 of the system, were one to exist; we point out that in [2] an 
asymptotic formula is given in the limit of a = oo and f3 large. 

In further detail, the motivation for the current paper, which completes the analysis of [4] and 
also concentrates on numerics related to the ILFD problem, is as follows: 

1. To compute the coefficients of the powers of p explicitly, at least for the lowest perturbation 
orders, so as to give a quantitative expression for the width of the tongues, for more general 
perturbations than those considered in [4]. 

2. To investigate numerically how large p can be for the analysis, which is carried out under 
the assumption that /jCI, to break down. 

3. To compute numerically the Arnol'd tongue diagram in the (co,p) plane in the case that the 
periodic part of the perturbation contains only one frequency, u. This allows us to obtain 
information for values of p where perturbation theory does not apply. On the other hand, for 
smaller values of p, the analytical results provide a check on the reliability of the numerical 
analysis. 

4. The same as 3, but in the case that the perturbation contains all integer multiples of w. it was 
argued in [4] that the width of all tongues would then be proportional to p and all the centre 
lines would be vertical. In particular we want to determine the constant of proportionality, 
i.e. the coefficient Aiu(p) in (1.2), and show that the higher the values of p and q in p = p/q, 
the lower the constant. 

The rest of the paper is organised as follows. In Section 2, we summarise definitions and lemmas 
from [4] which are needed in the remainder of the paper, and extend the analysis to more general 
analytical driving, possibly containing all harmonics. In Section 3 we concentrate on analytical 
results concerning the Arnol'd tongues, by gathering together all information which can be obtained 
to any order of perturbation theory. In Section 4, we describe the algorithms used to carry out 
the computations of the integrals appearing in the theory. In Section 5 we give and discuss the 
numerical results: after checking that they agree with the theory where the latter applies (small 
p) , we investigate how large p can be for the theoretical predictions to be reliably used. Finally in 
Section 6 we draw some conclusions. 

2 Preliminary analytical results 

We recall the results of [4], and extend them to more general perturbations. Numbered lemmas 
which we refer to in this paper are taken directly from [4], and all proofs are given there too. 
Reference to [4] is given only for proofs and technical details, the discussion below being quite 
self-contained. 

The system of ordinary differential equations that describes the ILFD can be put into the form 

u" + u' h(u) + k(u) + p^(u, u', t) = 0, (2.1) 
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with 

h{u) := 1 - 13 + 3/3u 2 , k{u) := u (a - /3 + (3u 2 ) , (2.2) 

and 

u', i) := w' (3u 2 - l) f(ujt) + u(u 2 - l) (f{ujt) + ujf{ujt)) , (2.3) 

where here and henceforth /' denotes the derivative of / with respect to its argument. The case 
f(t) = sint (and hence f'(t) = cost) was explicitly considered in [4]. More generally one can 
consider any analytic 27r-periodic function 

oo 

/(i) = ^jUn^, |/„|<*e-«M, (2.4) 

where the bound on the Fourier coefficients /„ — for suitable positive constants <!> and £ - 
follows from the analyticity assumption on /. For simplicity we confine ourselves to odd functions: 
considering functions whose Fourier expansion contains also cosines would overwhelm the analysis 
without shedding further light on the results. 

For /i = 0, (2.1) reduces to the Lienard equation [6, 10] 

u" + u' h(u) + k(u) = 0, (2.5) 

which we refer to as the 'unperturbed equation'. In order for it to have a globally-attracting limit 
cycle encircling the origin [10, 24] we require that a > (3 > 1 (this corresponds to the region 
of the parameter plane called design area in [7]). In that case, we designate uo(t) the solution 
to (2.5) corresponding to the limit cycle. Let T be the period of u (t) and let fi = 2n/T be the 
corresponding frequency: depends solely on a and (3. 

The unperturbed equation is autonomous, hence it clearly has the property that if uo(t) is a 
solution, then so is uo(t + T) for any constant T. Consequently, we can fix the origin of time so 
that uo(0) = Uq > and u' a (0) = 0. This has the effect of shifting the third argument of ^ by 
some time t n , so ^(u,u\t) becomes ^(it, u', t + t ) in (2.1). 

We also note that the symmetry properties of h{u) and k(u) guarantee that u n (t) has the 
property 

u Q (t + T Q /2) = -u (t) Vt £ R, (2.6) 

which in turn yields that the Fourier expansion of u n (t) contains only odd harmonics (lemma 2.1). 

It is convenient to rescale time by defining r = tot so that now has period 27T in its third 
argument. After rescaling, the differential equation becomes 

11- 
il+-u h(u) + -^k{u) + ^(u, ii,T + r ) = 0, (2.7) 

where a dot denotes differentiation with respect to r, tq = coto, and we have defined 

*(u, u, t) = \ [lou (3« 2 - 1) /(r) + u (u 2 - l) (/(r) + w /'(r))] . (2.8) 

id 

We have shown in [4] that if to is 'close' to p£lo/q, where e I are relatively prime, then 
the frequency f2 of the solution exactly equals qw/p. the system is said to be locked into the p: q 
resonance. How close ui has to be to pSlo I Q depends on \i and on the resonance itself — quantitative 
investigation of this 'closeness' is the aim of the present paper. 
Let p = p/q G Q. For u) close to p£l put 

11 °° 

- = —— +e(/i,T ), where e(/x, r ) = V ^ fe £ fe (r ). (2.9) 
lj pil £-< 
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Unlike [4], for the sake of convenience, here we make explicit the dependence of e on t . The 
perturbation calculation is then carried out by substituting the expression (2.9) for ui in (2.7) and 
expanding in powers of p. This results in 

oo 

H(u,ii, ii, p) := H (u, u,il) + p k H k (u,ii, t + r ) = 0, ( 2 -10) 

fe=i 

where 



iih(u) k(u) 



H (u,u,u) = u+ — ±± + (2.11a) 



g fc (M, it, r) = £ fc (r ) ( u h{u) + J + £fci(To)efe 2 (To) fc(u) 



ki+k 2 =k 



+ e fc _i(r ) 



ti (3u 2 - 1) ,f(r) + u (u 2 - 1) + /'(r)) 

+ X! £ fci( T o)£fc2(fo)w (w 2 - 1) /(r), fc>2, 

^1+^2 — ^—1 



(2.11c) 



where the last line of (2.11c) is missing for k = 2. 

In order to carry out the perturbation calculation to first order, we first write the unperturbed 
system in the form 

v h(u) k(u) . 
u = v, v = ^-^eG(v) (2.12) 

which has a unique 27r,o-periodic solution (u (t), vo(t)) such that «o(0) = 0. The Wronskian matrix 
of equation (2.12) is 

w{t)= (w 11 (t) Wl2 (r)\ (213) 
1 tun(r) I0i2 (r) ' 



and satisfies 



\w{t) = M(t) W(t), ( 1 \ 

\W(0) = 1, ^ ^-G( Uo (r),«o(r)) -G(u (r), ^(r)),/ 1 J 

Lemma 4.1 then states that a solution to equation (2.14) is obtained by setting 

rr e -F(r') 

wi2(t) := c 2 Uo(t), Wu(t) := CiU (t) At' , 2 (2.15) 
where J^(r) is given by 

fto:=-4- f Ar'h(u {r')), (2.16) 

the constant f £ (0,7rp) is chosen so that ti>n(0) = 0, while the constants Ci and c 2 are such that 
W(0) = 1 — it is shown in [4] that this choice can always be made. 

By defining r\ := uo(0), as in [4], and substituting this into (2.12), we find that 

n __ M-j£i»® . ( , 17) 

p SJq 

By Remark 4.2 in [4], we have, additionally, that c\ = —r\ and c 2 = l/r\. 
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We further define /o by p^ofo — (h), the mean value of h (u (t)), so that /o = F (2tt p) / (2n p) , 
and we write F(t) = f r + F(t), where F(t) is a 27rp-periodic function with zero mean. By 
lemma 1.2 one has /o > (cf. also [6]). 

Lemma 4.4 then states that there exist two 27rp-periodic functions a(r) and 6(r) such that 



Wu(t) = a(r) + e taT b{r), w 12 {t) = ca(r), 



(2.18) 



for a suitable constant c. In order to develop perturbation theory for a 27rp-periodic solution, with 
pel, which continues the unperturbed solution when p ^ 0, one writes 



u(r) = u (t) + ^2 V ku k{T), 



(2.19) 



fe=i 



where uo( T ) has period 2irp (and hence frequency 1/p). We have shown in [4] that there exist 
27rp-periodic functions Ufe(r) such that the perturbation series (2.19) converges for p small enough. 
The functions Ufc(r) are recursively defined (see equation (7.2) of [4]) as 

itfe(r) = iOn(r)ufc +w 12 (T)v k + / dr'c F(r,) [wi 2 (t)w) U (t') - Wu(t)w 12 (t')} *fe(r), (2.20) 

Jo 



with 



where 



*fe(r) := 



^ p k 'B w {u{t),u(t),t + r ) + G 2 (u(t),u(t)) 



L fe'=l 



(2.21) 



J fc 



G 2 («, u) := G(u, v) - G(u (t),v (t)) 

- (« - uo(r)) ^G(«„(r), t*(t)) - (v - «o(t))) ^G(uo(t), i*(t)) (2.22) 

and the notations means that we expand u{t) and w(r) according to (2.19) and, after taking the 
Taylor series of the functions Hk< , k' = 1, . . . , k, and G 2 , we keep the coefficients of all contributions 
proportional to p k . In (2.20), the initial conditions must be suitably fixed (again we refer to [4] 
for details), whereas can be set equal to zero (cf. remark 5.1 of [4]). 

Considering first order in p, we obtain the first order compatibility condition that has to be 
satisfied if U\{t) is to be periodic, i.e. (e F b^i) — 0, where ^i(t) = — i?i(u (r), v (t), t + r ). 
Expanding /(r) according to (2.4) and using (2.11b), this gives 



where 



and 



E>i 2v 



£i(ro) A + /„ L3,-i„ cos ^t + Sj2^ sin ^r ] = 0, 
i/=l 3 = 1 



1 /-27TP 

:= / dTC F(r) fe(r) 



u (r) /i(w (r)) + -rr-k(u {T)) 
pil 



1 K^t) . 

dr »„„ sint/r, 1 = 1,2 



o 



27rp 
1 

2?rp 7o 



2-rrp 



P 2 % 



B^iv '■= vpVl a B 22v , 
1;2, := — ^pilo^i!/, 



(2.23) 

(2.24) 

(2.25a) 
(2.25b) 
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with 

K^t) = c p ^b(r) P n vo(t) (3u 2 (r) - l) , K 2 (r) = c p ^b{r) u (r) (u 2 (r) - l) . (2.26) 
By setting D lu = - (B xlv + B 2Xv + B zlv ) and D 2lJ = - {B\2v + B 22v + B 32iy ), (2.23) then becomes 

£1 (r ) = — f v ( D\ v cos vtq + D 2v sin vtq) :=Di(t ). (2.27) 

V— 1 

By construction e\ has zero mean, so that either it is a non-constant function or it identically 
vanishes. For purposes of comparison with [4], in the following we shall shorten D\\ = D\ and 
D 21 = D 2l and also B t ji = Bij, which are the only relevant constants when / contains only the 
first harmonics v = 1 in (2.4). 

It is shown in Appendix B of [4] that A = —ripQo; hence, from (2.17), 

A - Wo ' ( 8) 

which provides an obvious means to check the numerics — by calculating A from (2.24) and 
comparing with (2.28). 

In [4] it is also shown how to go to higher orders; to any order k > 1 one finds the compatibility 
condition (e F b^h) = 0> where the function ^(r) is given by (2.21). 

The compatibility condition leads to all orders to equations like (2.27), which now read 

£fc(r ) =S> fc (T ), k>l, (2.29) 

for suitable functions Dk — strictly speaking in [4] only the case f(t) = smt is explicitly discussed, 
but one can easily work out the general case of / an arbitrary analytic function by following the 
same strategy. Note that, with respect to [4], here we have included the factor 1/A in the definition 

Of £>fc(T ). 

The width of the plateau corresponding to a given p (i.e. to a given resonance p : q such that 
p = p/q) can then be expressed as follows. First one defines 

oo 

T)(t ,p) =y"e fe Dfe(To), e max (p) := max S(<r ,p), e min (p) := min D(r ,p). (2.30) 

^— ' 0<T <27r 0<T <2-K 

k=l 

Then by setting 

Wmin(p) := TZ 7 T j ^max(p) := ~ pj 7 T ; (2.31) 

1 + pn £ max (p) 1 + P^O £min (P) 

the plateau corresponding to p is given by 

A/\ / \ / \ P ^0 ( £ max(p) — £min(p)) /- on \ 

Aw{p) := w max (p) - w min (p) = — — — — —77- 7-vT- ( 2 - 32 ) 

(1 + pil £min(p))(l + /O"0£max(p)) 

In other words, for u> G [w m i n (p), w max (p)], one has locking u> = pCl, if Q denotes the frequency of 
the output signal. For each such value of ui the initial phase To gets fixed to a value Tq such that 
1/u) = 1/pfio + ^(fi, Tq), according to (2.9). 

When the function £i(t ) in (2.27) does not vanish, then, if one further assumes that the 
second derivative of Si is non-zero at the stationary points (where the maximum and minimum 
are attained), the first order approximation is adequate. In other words, in such a case one can 
approximate 

£ max (/o)=A* max 2>i(t ) + 0(p 2 ), £min(p)=/" min £>i(r ) + 0(m 2 ), (2.33) 

0<T <27T 0<T <27r 
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and hence 

w m in(p) = P^o I 1 - P^oP max Si (r ) I + 0(p 2 ), (2.34a) 

y 0<r <27r y 

Wmax(p) = ptto ( 1 - P^oP min Si(t ) ) + C>(p 2 ), (2.34b) 
y o<to<27t y 

which gives a plateau of width 

Auj(p) ^ pA 1 uj(p) + 0(p 2 ), Aiw(p) := /r 2 ^ 2 , ( max £>i(t ) - min Di(t )) (2.35) 

\.0<To<27r 0<ro<27r y 

centred 'around' the value ui c (p) = pf2 - Since the function £i(to) has zero mean, this means that 
the corresponding Arnold tongue in the (w,p) plane emanates from the point w c (p) of the w-axis 
as a cone with axis along the vertical and angle 9{p) = #i(p) + 9 2 {p) such that 

tan 9i(p) = -p 2 nl min £>i(t ), tan6» 2 (p) = p 2 f^ max Di(t ). (2.36) 

0<ro<2-rr 0<ro<2-7r 

If / contains only one harmonic, say /„ = for \u\ > 1 in (2.4), then 9\{p) = ^(p), and max D\{tq) 
= A~ X ^J D\ + D 2 . Note that in such a case the second derivative of T>i equals ±f\/A when the 
first derivative vanishes. 



3 Arnol'd Tongues: analytical results 
3.1 First order contributions 

Let us consider the expression in (2.32) for the leading contribution to the width of the plateau 
when the first order contribution does not vanish. Then we neglect the high order terms and 
approximate 

Aw(p) « p,p 2 n 2 l Q(p), where Q(p) = max X>i(to) — min 2)i(to), (3-1) 

0<To<27T 0<To<27T 

Note that, to obtain Di(to) from (2.27), one must keep only the summands such that f v ^ 0. 

By writing Bjj„ according to (2.25), one uses that the Fourier expansions of the functions Ki 
contain only even harmonics (cf. Section 6 in [4]), i.e. 

Ki{r)= ^' T,P K IV ^Y.^ IJ ' T/Pk ^'y ( 3 - 2 ) 

v even 

Furthermore, as (2.26) shows, the functions Ki are analytic and hence the corresponding Fourier 
coefficients K iv i decay exponentially, i.e. for i = 1, 2 and for all i/'eZ one has \K iv *\ < Te - ^ 1 ^ 
for suitable positive constants T and £i. 

Hence by expanding Ki according to (3.2) and writing 

sini , r = J2 ^ iavT , cosvt= J2 \^ av \ (3-3) 

CT=±1 1 cr=±l 

one realises that one can have B^ v ^ only if there exist v' G Z such that K iv > y^ and 
v'q + avp = 0. If we assume that the first condition is satisfied for all even v' £ Z (numerical 
analysis ensures that such an assumption is reasonable — see figure 1), then the key condition is 
that there exist v'eZ such that 

2\v'\q = \v\ V (3.4) 
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with p, q relatively prime integers. When this happens one has 
1 



r U i/'eZ,(7=±l 
1v' -\-av p— 



^2 Ki(2 V ')Rja, i,j = 1,2, where R la = ^ R 2a = \, (3.5) 



and B 31v = vpVL^B^v, B?>iv = — vp£l®B 2 \ v . If the function / contains only the first harmonics 
(so that /„ 7^ only for \v\ = 1) then in (3.4) one has to consider only the case \u\ = 1. Thus, as 
discussed already in [4], one obtains q = 1 and p = 2\v'\, i.e. p must be even. This means that one 
finds plateaux of width O(p) only for resonances p : q with q = 1 and p e 2IN. 
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Figure 1: Fourier coefficients of the functions Jfi(r), x, and K'2(t), +, for a = 5, (3 = 4. The odd 
coefficients turn out to be zero (within the numerical error of ~ 1CP 11 ), according to (3.2), while all the 
even coefficients are non-zero and decay exponentially. The dotted lines are only to guide the eye. 

On the contrary, if the function / contains all the harmonics, the condition (3.4) has to be 
considered for all feZ, and one finds easily non- vanishing contributions to (3.5), e.g. by taking 
\v'\ = p and \v\ = 2q. Thus, for any resonance p:q one has a plateau which to first order is given 
by (3.1). From (2.23) one obtains 



Si (to) = 



1 



A P 2 n 2 



( K H"P) cosl ^ T o + K 2(u P ) siniAro) 



up even 



where we have defined 

K lv := J2 [ R ^ (^i 

<r=±l 

K 2v := ^ [ R 2c (Kl(-<Tv) + ^2(-<ti/)) - vpVLnRi <y k 2 (- C y V ) 

<T=±1 

Let us define v = min{^ > I : vp even} and v\ = minjV > uq : vp even}, and set 

2 \L\\k V0P {p)\. 



K v {p) := yl\K lv \* + \K 2v \*, Qo(p) = 



Then one obtains 



Q(p) = Qo(p) + o 



\A\ P *nl 

\.L\\K„ 0P (p)\ 



(3.6) 

(3.7a) 
(3.7b) 

(3.8) 
(3.9) 
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which inserted into (3.1) gives 



^ {p)K 2 J^L\f va \\ Kuap {p)l (3.10) 



where we have used that A = —n/pSlo, with the constant r\ independent of pfio- 

If one keeps the whole sum in (3.6) one finds, always in the first order approximation, 



|Au,(p)| < ?^ m ax^|/,||^)(p)| < 2j WT^ ^(l + ^O )|/,||%p)(p)|. (3.11) 

up even 



1^— 1 1 1 



Since p and </ are relatively prime the condition vp G 2Z can be satisfied only if \v\ > q and 
\ V P\ > P- Therefore for fixed p = p/q one has 

|Aw(p)| < /iCpVV^e-^, (3.12) 

where C is a constant independent of p. This shows that all the Arnold tongues have width 
proportional to p, but the constant of proportionality decays exponentially with p and q. Therefore, 
for fixed p, the union of all the Arnold tongues is O(p), and hence tends to zero when p — > 0, as 
expected from [11]. 

For instance, if f(t) = sini + ?7sin2t in (2.4), one has A(2n) = c(2n) p + 0(p 2 ) and A(2n — 1) = 
c(2n— l)r?/x + 0(p 2 ), for suitable constants c(n) independent of p and 7/. Therefore, for all integer 
resonances the plateaux are of the same order of magnitude — provided, of course, r) is large 
compared with p. 

3.2 Second order contributions 

When the first order dominates, the second order gives a correction which can be computed ex- 
plicitly. When the first order vanishes, the second order becomes the leading order (if it does not 
vanish too). 

To compute the second order one needs the function D 2 appearing in (2.29) for k = 2. The 
analysis in [4] and (2.29) show that 

{e p b^ 2 ) =Ae 2 (To) + (c p bZ 2 (-;T )) => D 2 (r ) = -j(e p bZ 2 (-; t )>, (3.13) 

where, by (2.11) and (2.21) with k = 2, one can write 

S 2 (t;t ) = S 2 (t;t ) + S 2 (t;t ), (3.14a) 
S 2 (t; t ) - -s 2 (to) [(a - (i) u (t) + /3u 3 (t)] - £i(t )« (t) (3u§(t) - l) f(r + t ) 

- £i (t ) («g(r) - uo(r)) (J^f(T + + f(r + D)) . (3.14b) 

S 2 (t;t ) = -«i(t)- — (mo(t),uo(t),t + to) - Mi(t)^t -(w (t), v (t), t + t ) 
1 9 2 G 9 2 G 

with H\(u,v,t) and G(u,v) given in (2.11b) and (2.12), respectively (we have explicitly used that 
G(u,v) is linear in v). 

Thus, to compute (3.13) one first needs the first order solution (ux,v\), with v\{t) = Ui(t). 
For k = 1 equation (2.20) gives 

ui(t) = co(t) (Q 2 (t) - Q 2 (0) - Qi(0)) - c6(t)Q!(t), (3.15) 
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where the functions Q\ and Q 2 can be read from equations (5.3) and (5.4) of [4], which we rewrite 
here for convenience, 

f T dT'e F ( T \(T')fi(r') =e^Q 1 (r)-Q 1 (0), f dr'c^( T ')6(r')*i(r') - Q 2 (r)-Q 2 (0), (3.16) 
Jo Jo 

and we are using that Q := (e F 6^'i) = and u\ = — cQi(O). Note that the functions u\ and v\ 
depend also on To; more precisely, by construction one has 

Ul ( r ) =EE e ivT ^e^ T+T ^U lvvl , U lvvi oc f Vl ; (3.17) 



v odd 



as easily follows by reasoning as in the proof of lemma 8.2 in [4], the only difference being that / 
can contain all harmonics. 

In particular, when £\ vanishes identically, then S 2 also is identically zero and (3.13) reduces 

to 

D 2 (r ) = - A — I drc^)6(r) — { [(^(r) ^(r) + — (3^(r) - l) )/(r + 

+ (3u§(t) - 1) f'(r + ro)] ui(t) + (3u§(t) - l) /(r + r ) «i(t) (3.18) 

+ ^(«b(r)/i"(«o(r)) + fc "y )} ) "iW + h '( u o(r)) «i (r) m (r) } , 

where = 6(3u, h"(u) — 6/3, and k"(u) = 6/3u (here, as well as for /, the prime denotes 

derivative with respect to the argument). 

By using the expansion (3.17) for m, one finds that the function 

D 2 (r ) = J2 C " T0 ®2^ = 2) 2 ,o + Yl ciUT °®^ ( 3 - 19 ) 



i^/0 



can be written in the form 



1 f P 

®2(to) = —Y< E &Tj VT/p e i{vi+V2){T+To) k v , VlV2 , U^e^H/j^ (3.20) 

z/ even 

for suitable coefficients h v i^ lV2 . Then one sees that only the coefficients k^2v')v 1 v 2 with 

2W\q=\u 1+ u 2 \ Pl LL^O, (3-21) 

contribute to (3.20). The term with v\ + v 2 = v' = gives the mean D2.0 of £) 2 , and requires no 
condition on p. This explains why the boundaries of the locking region are either O(p) — when 
the first order dominates — or 0(p 2 ) — in all the other cases. However, the width of the plateau 
arises from the variations of £> 2 , hence it is related to the terms in (3.20) with v 7^ such that 
(3.21) is satisfied. If there are any of such terms, then the function D2 is not identically constant, 
and therefore, in such a case, one has 



Au(p) = /rA 2 w(p) + 0{pf), A 2 uj(p) := p z nU max 2j 2 (t ) - min £> 2 (r ) (3.22) 

\0<i"o<27r 0<t <2-k J 

which replaces (2.35) when the first order vanishes. For instance if / contains only the first 
harmonics then (3.21) is satisfied for q = 1, p G IN and v\ = v 2 = ±1, which shows that the 
plateaux corresponding to odd p are of order p 2 — see [4] for further details. 
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3.3 Higher order contributions 



If one wants to determine the higher order contributions, the analysis above can be easily extended, 
even if it becomes much more complicated from the computational point of view. In general, if 
one writes 

oo 

Aw(p) = A fc w(p), (3.23) 

fe=i 

one finds 

k 

A k cj(p) = Y^ A fc w(p;n,...,iyfe), A k u(p;u 1 ,...,iy k )^e- 2 ^ p l[f 1/i , (3.24) 

vez i/i,...,i/|,ez i=i 

\vi + ...+v k \p=2vq 

so that, in order to single out the leading contribution to the width of the plateau, one has to 
compare the size of the perturbation parameter p with the amplitudes of the harmonics f v of the 
drive. 

Note that to all orders k the coefficients A k uj{p) decay exponentially in both p and q. Thus, 
every time the first order does not vanish it dominates, provided p is small enough. If on the 
contrary A k 'Lj(p) = for all 1 < k' < k and A k uj(p) =/= then one has Acj(p) — 0(p k ) for p small 
enough. 



4 Numerical computations 

4.1 Numerical solution of the ODE 

Since in general no closed-form solution to (2.7) with p — exists for (5 ^ 0, it is clear that this 
equation must be solved numerically. Furthermore, in order to approximate uq{t) and uq (t), the 
ODE must be solved for a sufficiently long time that the transient has, for practical purposes, 
decayed to zero. An effective initial procedure was found to be (i) solve the ODE from r = —t\ 
to t = 0, where t\ is large, using any standard method, for example, the Rungc-Kutta fourth 
order method; (ii) solve for a further small time T2 which is such that uo(t2) = and uo(i"2) > 0, 
again using the Runge-Kutta method, and additionally using bisection to find such that the 
first condition is met; (iii) solve from r = t-i to T3, where T3 is the smallest value of r which is 
greater than ti and for which, again, 110(13) = and uq(t^) > 0. Then an estimate of the period 
of uq(t), To, is T3 — T2 and an estimate of Uq is 110(72) ~ wo^)- 

In practice, these estimates can then be somewhat improved by solving the ODE assuming that 
a power series for Uq(t) exists, and computing this series around r — 0, using the initial conditions 
u (0) — [To, as estimated above, and «o(0) = 0. We can shift the origin of time from r 3 to zero 
because the ODE is autonomous. Typically, several power series need to be computed to cover 
the range r = to To, but the method has at least two advantages over Runge-Kutta. The first is 
that the error can be estimated by implementing a test on the coefficients of the power series, as 
set out in [5]; the second is that the Newton-Raphson method can be used directly on the power 
series for the solution around To to find the value of r for which uq(t) = 0, and hence to estimate 
T . The series used in practical computations had degree 30. 

Once accurate values of Uo and T have been computed, it is a simple matter to calculate a 
table of values of «o(r) and u (r) at r = ih, i = . . . M— 1 for some M G IN and for h = T /M > 
a given time-step. 

4.2 Interpolation 

A discussion of a suitable interpolation method is now in order. 
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In what follows, we will need to integrate functions of uo(t),Uo(t) and to do this we use an 
interpolation scheme from which such integrals can be computed directly. 

We start by discussing a scheme for interpolating from the values of uo (r) , iio (r) at discrete 
times ih, i = . . . M — 1, produced by the numerical ODE solver. 

Since u (t) is periodic, the interpolation scheme should reflect this — standard methods based 
on the Lagrange formula or Chebyshev polynomial interpolation are therefore not suitable. Instead, 
interpolation based on the function 

smKirr 

Ik{t) = — , 4.1 

K sin 7tt 

where K is an odd, positive integer, is used. This function is equivalent to a truncated Fourier 
series (see Appendix A) and has the properties that (i) Ik (t + 1) = Ik{t), so it is periodic (if K 
is even, the period is not 1 but 2); and (ii) 

/ 777 \ 

lim Ik(t) = 1 and Ik I T7 I = 0, m G Z, m not a multiple of K. 

T 7 > 5 \KJ 

Now let ie(t) = x(t + T ) be a periodic function of r with period T and set Xj = x(jT /K) for 
j = . . . K — 1. Then, defining 

K-l 

x{t):=Y,^Ik{t/T q -3/K), (4.2) 

we have, in the light of (i) and (ii) above, that x(kTo/K) = Xk = x(kTo/K) for k E Z. Hence, 
x(t) can be used to interpolate x(t) given the values of x(t) on a uniformly-spaced discrete set 
of values of r. In Appendix A we show that the error in the interpolation scheme described is 
O (e~ C2K ), for some positive constant C^. 

In practice, for r close to an integer, Ik(j) is best computed from a series expansion. Letting 
S = t — [t], with [t] being the nearest integer to r, we then use 



I K {8) = l-\{K*-l) 



{-Kdf ±(3K 2 7)M) 4 + ^(3^ 4 - 18if 2 + 31)M) 6 



+0(5 S ) (4.3) 



whenever |J| < ej. Since the computations are carried out to approximately 16 significant figures, 
we allow a margin of safety by choosing £j = 1CP 4 . 

The use of Ik(t) for interpolation has other advantages, amongst them that x(t) can be 
integrated in closed form, and so, by implication, the integral of x(t) for all r can be approximated. 
By defining 



r' K (T) = f 

Jo 



T , sin K-kt 
dr — ■ 



Sin 7TT 

it is easy to show that 



J'k+2(T) = J' K (T) + 2 f T dr cos(K + l)nr = J' K {T) + — j- sin(K + 1)ttT. 



Now, since K > is odd and J^T 1 ) = T, we have 

(if-l)/2 



sin 2i7rT. 



7T 

1=1 



Defining now J K (T) := At I k {t) = J' K (T)/K, we have 



T 1 (K - 1)/2 i 



1=1 
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Next define X(T) = f Q T Atx(t). Integrating (4.2) term-by-term, we obtain 

1 I \ / I 



i=0 



XiT)=T ^ Xi \j K (---) + J K [-)\, (4.5) 



where we have used the fact that Jk(t) is an odd function of r. In what follows, we therefore use 
X(T) to approximate jl' dr x(t). 

In a similar manner, defining Ek{Q,T) := J Q T dr e~^ T Ik{t), for constant (, it can be shown 
that 

T? (r T\ 1 - C ' CT , 2 (K \^ /2 C + c' CT (2»tt sin 2ittT - ( cos 2inT) 

Ek ^T) = ^^ + - ^ ^¥ ' (4 ' 6) 

Hence, X e «,T) := dre-^x(r) is given by 

X e ((, T) = To £ x% e-W* (| - 1) - £ K (-1)}. (4.7) 

4.3 Calculation of a(r), b(r) 

Before we can compute i»n(r), we need to find the unperturbed limit cycle, its period, T , the 
periodic function F(t) and the mean of F{t), /o- These are all straightforward: we solve the 
unperturbed equation (2.5) numerically as described in Section 4.1, obtaining C/ , To and the 
solution and its derivative over one period. Since u (t) is periodic, so is Ii(u (t)), and so we can 
use equation (4.5) to estimate F(t) for any r. From F(t) we can then obtain f , and hence F(t). 

Computation of wii(t) can now be carried out from equation (2.15), but is complicated by 
the fact that, for r = iT /2, ieZ, the integrand is singular and numerical integration techniques 
will break down. Singularity in the integrand, which is periodic, also prevents us from using 
equation (4.5). To discuss this further, let us define two subsets of R as S = U ie zSi with Si — 
[iT /2 — r c ,iT /2 + r c ], where r c T /2 is small and will be defined later; and I — R \ S. We 
will then use a power series representation for wh(t), t £ S, where the power series converges 
'usefully' (the error term is less than the maximum acceptable error) for |r| < r c , with Romberg 
integration [19], a standard numerical integration technique, being used for r £ I. 

It should be pointed out here that we do not need to compute f explicitly. Instead, we can set 
the lower limit of the integral to any convenient value, t; say, provided we add a suitable multiple 
of Uo(t); so our definition becomes 



wu(t) = M (r)fc 2 + u (t) / dr'ci , 2 , (4.8) 

Jn u o\ t I 



where the constant k 2 , which depends on r ; , will be chosen to ensure continuity. 

In practice, we only need to know u>ii(t) over a length of time consisting of two periods of 
«o(t), and so we calculate it for r £ [0, 2T ]: from Appendix A in [4], we know that wh(t) is 
well-defined even at t = 0, which we take to be our value of 77. 

We derive the formal power series for wh(t) by using the method of Frobenius to solve the 
differential equation (2.5), with initial conditions chosen so as to ensure that the solution is on the 
limit cycle. Thus, uo(0) = Uq, uq(0) = 0, from which we obtain the power series in r for uo(r) and 
hence, using term-by-term differentiation, for uq(t). The latter is 



Mr) = U (a - (3 + /3C/ 2 ) 



+ (1 - /3 + 3/3t/ 2 ) — 



+ 0(r 3 ) 
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By Remark 1 in [4], c\ = — r\, which is the coefficient of t in the above series for uq(t), and 
so ci = Uo(a — (3 + 0Uq). Using the series for uo(t) = Uq + J*J" (It'uq{t'), and term-by-term 
integration, we can also find power series for F(t), e~ F ^ and hence for the integrand e - ^ 1 "-* /uq(t). 
Integrating this series from to r term- by-term, and multiplying by the series for ciUq(t), we obtain 
Wu(t) =l-(l-0 + 3/3£/ 2 ) t/2 + (1 - 2a + /3 2 (1 - 3[/ 2 ) 2 ) t 2 /4 + 0(r 3 ). Finally, we apply the 
remaining condition on ibn, that is, wn(0) = 0, which forces the choice of k 2 in equation (4.8) to 
be such that -k 2 U (a - f3 + (3U§) - (l - /3 + 3/3J7 2 ) /2 = 0. This gives 

M 

<i(r) « 1 + £ R ^ + ( tM+1 ) > ( 4 - 9 ) 

where 2i? 2 = a - /3 + 3/3U 2 ,, 6R 3 = a - j3{l + a - 0) + 3/3(1 + 3a - 4/3)ug -I- 15/3 2 u^ and so on. 
This is a truncation of the series actually used for r e sq. Using computer algebra, it is possible to 
extend this series to at least the term of order r 10 , expressing each coefficient of r as a polynomial 
in a, (3 and Uo — that is, without assuming numerical values for these parameters — although the 
higher order coefficients become quite complicated. 

The series (4.9) can be used to estimate wh(t) for t S Sj, j > 0, provided (i) the value given by 
the series is multiplied by (— iye~jfo T o/2 an d (ji) j s chosen so as to ensure continuity across the 
boundary of Sj. The term (— 1)° in the correction factor arises as a result of the property of uq(t) 
in equation (2.6), and the exponential factor comes from the definition of wh(t), equation (2.15). 
Hence, w\\{t) is estimated as 

wu(t) = k 2 ii {r) + {-l) j e- jhTo/2 w^{T - jT /2) + O (r M+1 ) (4.10) 

for r £ Sj. 

A ; — ; " B 



I I + 



D 



-1 S H S H 

I 1 I 

a- Sr 



-1 -1 

I 

a- 



Figure 2: The cases A-E to be considered when calculating Wii (r). The subsets Sj = [jTo/2 — r c ,jTo/2 + r c ] 
and Sj+i are shown as thick line segments. It is assumed that wii is being computed at equally spaced 
time steps of width h, that wn(r — h) has just been calculated, and that ton(r) is now to be found. How 
this is done depends on the relationship of r to r — h — see text. 



In order to compute 6(r), we need to know wu(t) for t G [0, 2T ] — see equation (4.11) — and 
hence we calculate mjh(t) at equally spaced points 0, h, 2h, . . . , 2Xb, in that order, where 2To/h is 
an integer. The point r = is in so and so the truncated series is used here (with k 2 = 0). For 
t > 0, various different cases exist, and these are illustrated in figure 2, in which r is the time at 
which u>n is to be calculated, and it is assumed that it has already been calculated at r — h. 

• In case A, r € Sj, so the series is used, with the current value of k 2 . 
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• In case B, first wh(t + )/uo(t + ), where r + = jT /2 + r c , is calculated from the truncated 
series. To this is added a numerical estimate of dr'cie~ F ( T V"o( T ')' and the result 
multiplied by uq{t) to obtain an estimate of !»u(r). 

• In case C, roughly the reverse happens. Define t~ = jTq/2 — r c . Then numerical integration 
is used to estimate Wh(t~), from which k 2 can be found, by assuming continuity across the 
boundary r = t~ . Since the appropriate value of k 2 is now known, the truncated series can 
be used to estimate w\\(r). 

• In case D, compute as in C, followed by B. 

• In case E, straightforward numerical integration alone is used. 



In this way, wn(ih) is computed for i = 0, 1, 2, ... , 2T /h, and it is now a simple matter to 
extract a(r) and b(r) at the points r = ih, i = 0, 1, . . . , To/ft., so that their values at any r can 
be found by interpolation. From equation (2.18), we have that wh(t) = a(r) + c _ ^° r &(r) and 
a(r) = ju (t). Since a(r) and 6(r) both have period T , we have 

6(T)=e^ Wll(T) - Wll /l + ro) (4.11) 

and, knowing Wh(t) for r £ [0, 2T ], we can now compute b(r) for r e [0, To]. Having computed 
6(r), we can use for instance the method of least squares to estimate the value of 7: that is, we 
find the value of 7,7, that minimises 

T / h 

[wu(ih) - e~ folh b{ih) - jii (ih)] 2 , 

i=0 



which is 

D6o(i/i) [w n (ih) - e-f° ih b(ih)] 
where the sums go from i = to To/h. This completes the calculation of a(r) and 6(r). 



4.4 Illustrative results 

Illustrative results are now given for the case a = 5, = 4 and /(r) = sinr. All the computations 
were carried out using double precision arithmetic. For interpolation, K — 151 equally spaced 
points were used; the series for Ik (8) was used if \S\ < £j = 10~ 4 . In series (4.9), M — 10. In the 
definition of sj, r c = 10 -2 , and the fractional accuracy chosen for Romberg integration was 10~ 12 . 
With these parameters, and p = 2, we find T w 3.698939867513906, U Q « 0.979106186033891, f w 
0.757499334158 and 7 = -54.855909271256. Having computed o(r) and &(r), we can then estimate 
A, first of all from equation (2.24), using Romberg integration: this gives A — 16.0813516305191. 
Using equation (4.5) to carry out the integration, we obtain A — 16.0813516305189. Furthermore, 
we have from equation (2.28) that A = -r lP n a = U a (a - + 0U%) /{p9. ) = 16.0813516307791. 
These estimates agree with each other to 11 significant figures, thereby verifying the numerical 
techniques used to obtain them. 

The calculation of B\\ . . . B32 and D\, D2, defined after equation (2.27), now follows straight- 
forwardly from equations (2.25). The only point to note is that these integrals can be zero, which 
gives problems in the error control scheme used for numerical integration. To overcome this, the 
integration is done in two parts, from to 2npz and from 2irpz to 2irp, with z approximately, but 
not exactly, one half. We then find, for the above parameters, that, when p = 1, D\,Di w 10~ 12 . 
On the other hand, with p = 2, wc find D 1 w 8.11989 x 10~ 2 and D 2 « -5.20174 x 10 _1 ; for p = 4, 
D 1 = -3.79022 x 10~ 2 , D 2 = 2.74434 x 10" 1 . 
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5 Numerical results 



To extend our analytical results to large values of p we compute a set of Arnol'd tongues numerically. 
We make several comparisons between the theoretical predictions and the computational results, 
which provide information about the computational accuracy and the range of validity of some 
theoretical estimates. 

5.1 Arnol'd tongues 

We computed the Arnol'd tongues of system (2.1) for the 15 strongest resonances using the algo- 
rithms from [21] with two types of forcing: Harmonic forcing with 



containing all harmonics. Note that / in (5.2) is smooth and /(r) e [—1, 1]; hence, it can be used 
as a direct replacement for the sine function in (5.1). The relative strength of the harmonics can 
be adjusted by varying A, since one has /„ = <j>(A) A~^l, with $(A) = (A 2 - l)/2iA. 

The results of both computations for the parameter values a = 5 and j3 = 4 are illustrated 
and compared in figure 3. Note that in this case we have £1 = 1.698645 . . . and, hence, u) c (2) = 
3.397290 . . . and w c (4) = 6.794580 ... As explained in detail in [21] these tongues are computed by 
continuation of so-called constant-/^ cross sections, starting at the tips. To facilitate our subsequent 
computations of the order of contact and opening angles we started with an extremely small 
continuation step size to obtain a large number of points very close to the tips for later fitting. For 
each tongue we performed 150 continuation steps. The computation of most tongues terminated by 
either reaching the computational boundary p = 3.5 or by exceeding the maximal number of 150 
continuation steps. However, the computation of some tongues, most notably of the 2:1 tongue, 
seems to end due to limitations of the algorithm we use; see [21] for more details. We did not 
pursue a further investigation, because we are mainly interested in the size and location of the 2:1 
and 4:1 tongues for moderate p and in investigating the behaviour at the tips of all tongues, for 
which we obtained sufficient data. 

The bifurcation diagrams shown in figure 3 are clearly dominated by the strongest resonances 
occurring for p = 2 and p = 4. A continuation of the frequency-locked sub-harmonic solutions 
along the centre-lines lu/Q,q = 2 and u>/f2o = 4 inside these tongues revealed that these solutions 
remain attracting for p < 3.5 at least along these centre-lines. On the other hand, we observe at 
the left-hand boundary of the 2:1 tongue that this tongue overlaps with other tongues. Hence, in 
these overlapping regions we might find multi-stability. To the right-hand side of the 2:1 tongue 
no such phenomenon is apparent in these figures. 

The small plots to the right of figures 3 (a) and (b) show enlargements of the tip of the 3:1 
tongue illustrating the effect of forcing with and without all harmonics present as predicted in 
Section 3.1. In figure 3 (a) we observe a high-order (quadratic) contact of the boundaries of the 
tongue, while in figure 3 (b) the two boundaries intersect transversally. Note that the slight shift to 
the left of u>/Qo = 3 is due to the discretisation error of the periodic solutions; see [21] for technical 
details. It is remarkable, however, that our computations accurately capture the predicted high- 
order behaviour despite this approximation error, which is orders of magnitudes larger than the 
width of most tongues at their tips. 

To quantify our findings and for comparison with the analytical predictions we developed a 
simple adaptive non- linear fitting algorithm for the width-function Aui(p) to a monomial Aw(p) = 



f(t) = sin t 



(5-1) 



as considered in [4], and the forcing function 




A > 1 



(5.2) 
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(a) 



0.5 





co/Q 

Figure 3: Some Arnol'd tongues in the (uj/Qo, p)-plane for the ILFD for a = 5 and f) = 4 with (a) 
/(t) = sinr and (b) /(r) given by (5.2) with A = 2. For these parameter values we have O = 1.69864489, 
To = 2%/Q.q = 3.69893987, Ap = 2.78668166, Di = 0.00703534 and D 2 = -0.0450695. 



1.5 
1 

0.5 




Figure 4: Comparison of linear (dashed) and non-linear (solid) fits for a simple test example. 



a/j, with a and b unknown. One might argue that one could use a linear fit to logarithmic data of 
the form ln(Aw(p)) = In a + bln/i to compute estimates for a and b. However, this leads to biased 
estimates as we illustrate in figure 4, where we compare the results of a non-linear fit (solid) with a 
linear fit (dashed) to the function y — a exp(bx). Only the non-linear fit is a useful fit to the data as 
figure 4 (a) clearly illustrates, the linear fit here being biased towards lower function values. There 
are several reasons why a linear fit to logarithmic data is inappropriate, the most important ones 
being that the two least-squares residual functions \\Y — aexp(6AT)||| and || lnF — (lna + MQHl 
have different minimisers, and that Y and In Y do not have the same error distribution. Another 
suggestion could be to compute a linear fit to a polynomial p n (n) = a + a\\i + • • • + a n \i n of 
sufficiently high order n. However, we found that this leads to a least-squares problem that is so 
ill-conditioned that round-off errors become amplified to order one, that is, the fitted coefficients 
are essentially meaningless. A way out is to use orthonormal polynomials as base functions instead 
of monomials of the form au\i k '■ However, since our non-linear fit worked sufficiently well we did 
not pursue this further. 
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p:q 


/(t) = sin(r) 


/( T ) A^ + l-2Acosr' A 2 


Au) (p/q) /K fit iVfit 


Aui(p/q) /K fit iVfit 


1:4 
l.o 

2:5 
1:2 
3:5 


2.770 x 10" 7 At 8 076 8.763 x KT 1 48 
i v in - ^ .,6.052 Q71 w in -1 /in 

7.554 x 10~ 5 /i 5 038 8.294 x KT 1 36 
5.188 x 1(T 4 a* 4 ' 014 6.529 x KT 1 40 
3.093 x 10" 8 a* 10 19 8.839 x 10 _1 14 


5.571 x 10" 4 /i 1 001 4.481 x 10 -3 72 
z.y/i x iu (i o.yoz x iu to 
7.128 x 10~ 3 a* 1001 7.913 x 10~ 3 73 
1.786 x 10~ 2 /i 1 001 1-068 x 10 -2 74 
1.057 x 10" 5 Af 1001 7.612 x 10~ 4 58 


2:3 

1:1 

4:3 
3:2 


4.096 x 10" 3 fj 3 - 003 2.662 x 10 _1 49 

9 4DQ y 1 D~ 6 ,,8.016 7 4^9 v i ri — l 1 f; 

Z.IUy A ±U fJ/ 1 A J.U 111 

4.828 x 10~ 2 /i 2 000 9.099 x 10~ 2 76 
9.979 x 10~ 3 n 3 001 1.540 x 10 _1 47 
1.112 x 10" 3 /x 4 014 5.338 x 10 _1 38 


4.780 x 10" 2 /i 1001 1-612 x 10~ 2 74 
c ini y in~ 5 ,,o.999o 7oir v in~ 4 

O.lUl A ±U fJ/ 1 .U1U A ±U OO 

1.467 x 10" 1 /i 1 ' 001 3.179 x 10~ 2 76 
5.016 x 10~ 2 At 1 001 1-097 x 10~ 2 70 
1.688 x 10" 3 /i 1 001 1-973 x 10~ 3 61 


5:3 
2:1 
5:2 
3:1 
4:1 


5.049 x 10" 5 A* 5 971 2.588 x 10 _1 16 
7.556 x 10 _1 /i 1 000 9.684 x 10~ 2 80 
1.595 x 10~ 2 /i 3 971 2.039 x 10" 1 35 
1.024 x 10" 1 /i 1 " 9 8.960 x 10~ 2 77 
7.957 x 10 _1 At 09999 9.719 x 10 -2 80 


2.854 x 10" 5 At " 90 5.950 x 10 ~ 4 55 
6.280 x 10~ 4 /i 1 000 9.024 x 10 ~ 2 79 
1.832 x 10~ 4 At " 90 4.539 x 10~ 4 50 
1.331 x 10~ 2 /i ' 9994 2.4 10 x 10 ~ 2 73 
5.968 x 10 _1 /i "" 9.757 x 10~ 2 79 



Table 1: Leading contributions to the plateau widths corresponding to the main resonances as they appear 
in figure 3 from left to right. These coefficients were obtained by fitting the monomial ap b to the numerically 
computed values for Auj(p/q) over the interval fi <= [0, //fit] on Nfn data points. 



For our computations we used the weighted least-squares residual function 

F(a,b) := \\W(p) (Au(p) - ap b 



MW 2 



where we used the weight-function W(p) = 1/p to penalise errors closer to the tip a* = 0. The 
numerical data is normalised to the unit square, that is, we compute a fit to Aw(p)/ max{Au;(p)} 
and pj max{/i}. We computed estimates for a and b by applying Newton's method to the equations 
dF(a,b)/d(a,b) — and rescaled the computed coefficients to fit the original data. The size of 
the fitting interval a* S [0, pm] was computed adaptively. We started with an initial fitting interval 
/ifit = max{0.1, argmax{Aa;(A>) < 10~ 4 }}, that is, we used either pfn = 0.1 or the largest value of p 
such that the width of the tongue was less than 10~ 4 . The fitting interval was accepted if the least 
squares error F(a, b) for the normalised data was less than 1CT 3 and reduced successively if the 
error was larger. Within the fitting interval we excluded points for which Au>(p) was zero within 
numerical accuracy. Table 1 summarises our results for both forms of forcing. Each row states 
the fitted monomial representing the leading order term together with the fitting interval and the 
number Nfn of data points this monomial was fitted to. These computations agree extremely well 
with the theoretical predictions and also verify that our fitting algorithm is suitable to capture the 
leading-order behaviour accurately. 

First of all, our computations of the width of the locking intervals are in alignment with our 
theoretical results 

A ^ = { 0{fi) for p Odd!' 
for harmonic forcing (5.1), and with 

Aw(p/«) = O(p) 

for general forcing (5.2) containing all harmonics. Furthermore, for the main Arnol'd tongues 
corresponding to the resonances 2:1 and 4:1, equation (2.27) valid for harmonic forcing reduces to 

2>i(t ) = -^i D i cost + L>2 sin t ) , (5.3) 
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where, for the 2:1 resonance, A = 16.0814, D 1 = 8.11989 x 10~ 2 , and D 2 = -5.20174 x 10" 1 ; see 
also Section 4.4. An easy computation gives 



max Di(t ) = - min Si (to) = M, M = 0.0327381. (5.4) 

0<ro<27r 0<to<2tt 

For example, for the boundaries of the tongue corresponding to the 2:1 resonance we find from 
table 1 that (tan#i(2) + tan6» 2 (2))/2 w 0.7556/2 = 0.3778 in agreement with (2.32), which gives 
p 2 Q%M = 0.37785. Similarly, for the 4:1 resonance wc have A = 32.1627, D 1 = -3.79022 x 10~ 2 , 
D 2 = 2.74434 x 10 _1 , giving M = 8.6137 x 10~ 3 . Wc compute that p?Q%M = 0.39766, again in 
agreement with the result in table 1 that (tan#i(4) + tan<? 2 (4))/2 = 0.7959/2 = 0.39795. 

Also for the secondary resonances p: 1, with p odd, the agreement between the numerical results 
and the analytical predictions is satisfactory. The second order computation, performed according 
to the analysis in Section 3.2, gives, for the 1:1 and 3:1 resonances, the values A(l) w 4.8246 x 10~ 2 
and A(3) « 1.0269 x 10" 1 , to be compared with the values 4.828 x 10~ 2 and 1.024 x lO" 1 in table 1. 

For the all-harmonics forcing (5.2), the plateau widths as given in table 1 are consistent with 
the scaling law (3.10). If we fix p to be an odd integer, then v$p = v^pjq = 2p, hence v$ = 2q, and 
we find \f Vg \ \K Vg p(p)\ = $(A) \~ 2q \K 2 p(p/q)\- If, on the contrary, we fix p to be an even integer 
then v p = v^pjq = p, hence = q 1 so that we obtain \ f Vo \ \K Vop {p)\ = $(A) \~ q \K p (p/q)\. When 
inserted into (3.10), this leads to 

J cp/(q2^) for p odd, 
|Aw(p)| ~ <^ (5.5) 
I cp/(qz q ) for p even, 

with the constant c = c(p) independent of q. The constant c can be computed using the theory. 
A comparison between (5.5) and (3.10) gives c(p) = ln2rio\fi\~ 1 &(\) p\K U()P (p/q)\, with v$p = 2q 
for odd p and v^p = q for even p. For p = 1, 2, 3 we compute c(jp) = 0.82, 1.64, 0.11, respectively. 
These estimates are consistent with our numerical data in table 1. Fitting our data to function (5.5) 
we obtain the numerical estimates c(p) = 0.5867, 1.255, 0.05326, respectively, which is in good 
agreement considering the limited numerical accuracy and that (5.5) is valid for q — > oo. Figure 5 
shows a comparison between the theoretically and numerically obtained width functions. 




Figure 5: Tongue widths |Au(p)| for fixed p and varying q. The black curves are plots of (5.5) using the 
constants c from the numerical data and the grey curves are plots of (5.5) using the constants c predicted 
by the theory. The discrepancies are due to the small values of q are are expected to be asymptotically 
zero for large q. 
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5.2 Width of plateaux as a function of ct and (3 



Of practical importance is the rate at which the width of a given locking intereval increases with 
increasing p. The locking regions, i.e. the Arnol'd tongues, are cone-shaped and the vertical angle 
of the cone, 29\ (p) , which is a measure of width growth rate, depends on the parameters a and 
0. This angle can be computed from equation (2.36) for a given p, and typical results are given in 
figure 6 for p = 2 and a variety of values of 0, with a G (0, 10]. 




i i i i i i i i i i i 

0.0 2.0 4.0 6.0 8.0 10.0 

a 



Figure 6: A plot of 2#i(2), the opening angle (degrees) of the Arnol'd tongue for p = 2, against a, for 
various values of 0, with a > > 1. The angles have been computed numerically from equation (2.36). 



6 Conclusions 

In this paper we have investigated both analytically and numerically the structure of the Arnol'd 
tongues for a resonant injection- locked frequency divider (ILFD). It is the natural extension of 
the analysis performed in [4] , where we analytically proved the experimental and numerical results 
contained in [15, 16] by providing explicit formulae for the width of the plateaux appearing in 
the devil's staircase. More precisely in [4] we found the following result. Denote by ui and Vt 
the frequencies of the driving signal and of the output signal of the ILFD, respectively, with p 
the driving amplitude. Then, if for p G Q we call Au>(p) = {uj : ui/Q = p} the width of the 
corresponding locking interval, we showed that Aui(p) satisfies Aw(2n/fc) = 0(p k ) and Acu((2n — 
l)/fc) = 0(p 2k ) for all k,n G IN such that 2n/k and (2n — l)/k are irreducible fractions. In 
particular this implies that the largest plateaux correspond to even integer values of the ratio w/fi. 

In this paper we have extended the above results: we studied the system of ordinary differential 
equations (2.1), (2.2), (2.3), which describe the ILFD, with a more general driving term in the form 
of any analytic periodic function (we confined ourselves to functions containing odd-harmonics only 
in order to make the analysis more transparent and yet without any significant loss of generality). 
In [4] we used f(t) = sin(t) (one harmonic only), as in [15, 16]. Here, we studied the locking 
intervals Aco(p) by using in (2.3) a 27r-periodic function of the form f(t) = f^smi/t, with 

\f u \ < <I> c ^ I ^ I (by analyticity) . We found that, for any p = p/q G Q, with p, q relatively prime 
integers, the key condition for the existence of the locking u> = pfl (and hence of a plateau), is 
that there exists v such that /„/0 and 2\v'\q = \i/\p and some v' G Z. This condition is certainly 
satisfied if, for instance, \v'\ — p and \v\ — 2q, provided f 2q ^ 0, or = 2p and \u\ — q, provided 
fq ^ 0. Thus, for any resonance p: q one has a plateau Auj(p/q) which to first order is given by 
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(3.1) and (3.6). In particular, to leading order, the width of the Arnold tongues is expressed as 
Ao>i(p) w 2/xpfio|fi| _1 |/;y | \K Uop (p)\, where p = p/q, fig and fi are constants depending on the 
unperturbed system (but not on the driving) and v a > 1 denotes the integer which provides the 
leading coefficient in the sum (3.6). Note that the formula reduces to the one obtained in [4] — as 
it should — if f(t) = sin(i): in that case /„^0 only for v = 1, so that q = 1 and p 6 2IN. 

Moreover, by keeping in (3.6) the whole sum, we obtain |Ac^(p)| < pC p^^e"^ 3 ^-^ where 
C is a constant independent of p and q, thereby showing that the Arnold tongues have width 
proportional to p, but with proportionality constants which decay exponentially with p and q. 

We have also computed analytically the contribution of the second order, namely the coefficient 
of p 2 . In this case one needs to compute the first order solution (wi(t), tti(r)), with ui(t) given 
in (3.17), which rather complicates the analysis. We found Auj(p) = p 2 A 2 Lo(p) + 0(p 3 ), which 
replaces (3.1) when the first order vanishes. For instance, if / contains only the first harmonics 
then the condition for locking onto a p:q resonance becomes: 2\v'\q = \v\ + vi\p., with f Vl f V2 i= 0, 
is satisfied for some v' e Z. This shows that when f(t) = sint, as in [15, 16, 4], the plateaux 
corresponding to odd p are of order p 2 . 

Higher order contributions can in principle be computed with a very similar strategy (see 
(3.23) and (3.24)); the important point to notice is that to all order k the coefficients Akto(p) 
decay exponentially in both p and q. Naturally higher order terms become dominant when all the 
terms of smaller order vanish. 

To complete our investigation, we computed the functions ®i(tq) and 2)2 (to) numerically, 
from which the tongue widths Au(p) and Aw 2 (p) can be calculated, via equations (3.1) and (3.22) 
respectively. Some of the techniques required to carry out this computation are described in 
Section 3. We then computed a set of Arnold tongues, which was sufficiently large for testing 
the numerics on the basis of the theoretical predictions. In particular, we computed the width of 
the tongues for two types of forcing: (i) only one harmonic and (ii) all harmonics present in the 
Fourier expansion. Our computational results are in excellent alignment with the theory as stated 
above, which supports our belief that the locking charts in figure 3 are accurate. These two charts 
clearly demonstrate the dominance of the 2:1 and the 4:1 resonances. Furthermore, a comparison 
indicates that the location of the tongues is robust under generic perturbations; the differences in 
the shapes of the tongues are small. 



A Error in the interpolation scheme 

Starting from equation (4.1), we expand the sine functions in terms of complex exponentials to 
obtain 

M*) = ^ E (AT) 

\j\<(K-l)/2 

Setting t = t/T 0} so that the scheme can be used to interpolate a periodic function x(t) of arbitrary 
period T in terms of r, we have the Fourier expansion 

x(t) = J2 ®ne 2in7Tt . 

We interpolate x(t) by 

K-l 

x{t) - J2 x(j/K)I K {t-j/K) = Yl Pine 2 ™**, (A.2) 

3=0 \m\<(K-l)/2 

where the last inequality follows from equation (A.l). In order to determine how well x(t) is 
approximated by x(t), we need to compare a m with (3 m . Substituting for x(j / K) and lR(t) in 
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equation (A. 2) and rearranging, we find 



K-l 



x(t) = &M \ E X ^ e 2i ^- m ^/ K 

\m\<(K-l)/2 3=0 

where the term in braces is equal to (3 m . The sum over j is equal to [1 -c 2l ( rl_m ^]/[l _ e 2i(™-mW-^] 
provided that (n — m) ^ pK, peZ, and is equal to K otherwise. Hence, 

J_ e 2i(n-m)jir/K = J n - m ^ pK 

K ^ I 1 n — m = pK. 

3=0 ( - 

Hence 

Pm = ^ a m+pK = a m + a m - K + a m+K + ■ ■ ■ 



Now, since x{t) is the solution of an ODE with analytic coefficients, it is itself analytic, and so, for 
all n, \a n \ < Cie~ C2 ' n ', where C\, C2 are positive real constants. Thus, 

10m -ami < C\e^ K 

and hence, by choosing K sufficiently large, the interpolation error can be made as small as we 
please. 
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